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ABSTRACT 

With the diffusion approximation, we study the cascade and damping of 
Alfven-cyclotron ffuctuations in solar plasmas numerically. Motivated by wave- 
wave couplings and nonlinear effects, we test several forms of the diffusion tensor. 
For a general locally anisotropic and inhomogeneous diffusion tensor in the wave 
vector space, the turbulence spectrum in the inertial range can be fitted with 
power-laws with the power-law index varying with the wave propagation direc- 
tion. For several locally isotropic but inhomogeneous diffusion coefficients, the 
steady- state turbulence spectra are nearly isotropic in the absence of damping 
and can be fitted by a single power-law function. However, the energy flux is 
strongly polarized due to the inhomogeneity that leads to an anisotropic cas- 
cade. Including the anisotropic thermal damping, the turbulence spectrum cuts 
off at the wave numbers, where the damping rates become comparable to the 
cascade rates. The combined anisotropic effects of cascade and damping make 
this cutoff wave number dependent on the wave propagation direction, and the 
propagation direction integrated turbulence spectrum resembles a broken power- 
law, which cuts off at the maximum of the cutoff wave numbers or the ^He 
cyclotron frequency. Taking into account the Doppler effects, the model can nat- 
urally reproduce the broken power-law wave spectra observed in the solar wind 
and predicts that a higher break frequency is aways accompanied with a greater 
spectral index change that may be caused by the increase of the Alfven Mach 
number, the reciprocal of the plasma beta, and/or the angle between the solar 
wind velocity and the mean magnetic field. These predictions can be tested by 
future observations. 
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1. Introduction 

Turbulence is ubiquitous in the universe and plays important roles in our understanding 
of many natural phenomena (Kolmogorov 1941; Iroshnikov 1984; Kraichnan 1965). It occurs 
in highly non-equilibrium systems, where the microscopic viscous and/or resistive dissipation 
processes can not effectively convert the free energy into the internal energy of the fluid. 
Such systems usually have very high Reynold numbers and the free energy is stored in the 
large scale motion and/or magnetic fields. For plasmas, the free energy initially may also 
be stored in non-equilibrium distributions of charged particles. Turbulence is generated 
through a variety of instabilities related either to the large scale (magneto-)hydrodynamical 
processes or the microscopic collective plasma effects or plasma physics processes. These 
aspects have been extensively investigated with (magneto-)hydrodynamics and/or plasma 
physics theories. 

In astrophysics, most turbulence is carried by magnetized plasmas. It is responsible for 
distributing energies among different components of the plasmas, which may result in distinct 
emission characteristics or other observable features. Observations of these radiations can 
then be used to study the corresponding astrophysical sources (e.g. Liu et al. 2004; 2006 
on ^He rich solar flares). Many of the astrophysical plasmas we are interested in are made 
of charged particles. Depending on the turbulence energy loss rate to these background 
particles, its energy evolution is usually separated into two phases or ranges: cascade and 
damping. For the former, energy is usually transferred from large scales to small scales, and 
the turbulent motion has weak coupling with the background particles so that there are no 
significant energy exchanges between them. The corresponding spectral range is called the 
inertial range, where the energy flux is independent of the spatial scale. 

In the damping phase, there are strong couplings between the charged background par- 
ticles and turbulent motion, and the background particles are energized. For coUisional plas- 
mas, where the Coulomb collision timescales are much shorter than other relevant timescales, 
the background particles very quickly reach a thermal equilibrium. The energy partition in 
such systems is relatively simple, and its studies focus on other aspects, such as radiative 
processes, radiation transfer, ionization, and elements mixing. Much astrophysical turbu- 
lence, however, is carried by coUisionless plasmas, where the Coulomb collision timescales are 
long. Although it is generally accepted that the particle distributions in these plasmas are 
determined by the coupling of these charged particles with the turbulent electro-magnetic 
fluctuations, the details of these interactions are not well understood. Given the high degree 
of freedom to characteristize the particle distributions, it is usually assumed that the back- 
ground particles are in thermal equilibrium at least with their kinds. The damping therefore 
occurs on the smallest spacial scales and has been an essential part of plasma physics theories 
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for these kinds of coUisionless thermal plasmas (Andre 1985; Gary & Borovsky 2004). The 
damping of fast-mode waves by isotropic power-law electron and proton populations under 
typical solar-fiare conditions was recently studied by Petrosian et al. (2006). Large scale 
waves can be damped by high energy particles through resonant wave-particle couplings in 
the case. 

The cascade has been an essential element in all kinds of turbulence studies. The 
highly nonlinear nature of turbulence makes this a very challenging object. Nevertheless, the 
energy transfer of isotropic hydrodynamical turbulence, the simplest form of turbulence, has 
been described reasonably well with the Kolmogorov phenomenology that assumes a scale- 
independent self-similar cascade process. This assumption, in combination with the fact that 
turbulence energy is usually injected at large scales and dissipated at small scales, leads to 
the famous Kolmogorov power spectrum with a power-law spectral index of 5/3 in the inertial 
range (Kolmogorov 1941). Most astrophysical phasmas carry magnetic fields so that large 
scale magnetohydrodynamical (MHD) waves and small scale kinetic plasma fluctuatioins 
may be excited. The wave propagation effects on the cascade of MHD turbulence were first 
discussed by Iroshnikov (1963) and Kraichnan (1965). Although MHD effects are expected 
to introduce anisotropy to the system that can affect the couplings between the turbulence 
and background particles significantly, these preliminary investigations assumed that the 
turbulence is isotropic and reduced the cascade to a 1-D problem. The reduction of the triple 
correlation time to the wave period gives rise to a power- law spectral index of 3/2 (Zhou 
& Matthaeus 1990). These are called the Iroshnikov-Kraichnan (IK) phenomenology. Since 
the Kolmogorov cascade time scales as A;~^/^, where k is the wave number, (the IK cascade 
time scales as k~^/^) and the periods of MHD waves scale as k~^, the wave propagation 
effects become more important at smaller spatial scales and therefore play critical roles in 
the energization of background particles. 

To study the anisotropic effects induced by the presence of a large scale magnetic field 
and their implications on the plasma heating and particle acceleration, MHD turbulence in 
a uniform medium has been treated as an ensemble of linear wave modes (e.g. Petrosian & 
Liu 2004). Its nonlinear nature is revealed in the wave- wave couplings. Although this quasi- 
linear treatment may not be valid for strong turbulence, it is certainly a good approximation 
when the magnetic field fiuctuations are much smaller than the large scale field and the wave 
periods are much shorter than the eddy turnover timescales. Significant insight of properties 
of Alfven and magnetosonic turbulence has been obtained recently through this approach 
(Sridhar & Goldreich 1994; Goldreich & Sridhar 1995; Galtier et al. 2000; Chandran 2005). 
However, the component of Alfven turbulence excitations nearly perpendicular to a large 
scale magnetic field is dominated by nonlinear effects (Montgomery & Turner 1981). This 
indicates inherent limitations of treating MHD turbulence as a spectrum of waves. As shown 
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by Montgomery and Matthaeus (1995), linear MHD wave modes do not give a complete 
description of turbulence excitations and their coupling, and treating MHD turbulence as an 
ensemble of linear wave modes may miss some critical nonlinear effects (Ng & Bhattacharjee 
1996). 

Based on the dominance of nonhnear or wave-propagation effects, Oughton et al. (2006) 
recently separated the Alfven turbulence into two interacting parts: quasi-two-dimensional 
and wave-like fluctuations. The quasi- two- dimensional component characterizes the nearly 
perpendicular excitations and may be described with a quasi-2-D Kolmogorov phenomenol- 
ogy. The wave-like fluctuations may be described with an IK phenomenology with prop- 
agation direction dependent wave periods. The cascade of turbulence then also depends 
on the assumed couplings of these two components. It is interesting to note that the two 
components are separated by the critical balance between linear wave periods and nonlinear 
eddy turnover timescales, and the incompressible Alfven turbulence described by Goldreich 
& Sridhar (1995) appears to be appropriate for the quasi-two-dimensional component that 
is dominated by nonlinear effects. Due to the suppression of cascade by wave propagation 
effects, Alfven turbulence with wave periods proportional to the parallel component of wave 
vectors cascades preferentially in the direction perpendicular to the large scale magnetic 
field. The wave-like fluctuations need longer time to develop, which may explain the expo- 
nential cutoff of the power spectrum in the direction parallel to the large scale magnetic field 
observed in MHD simulations (Cho, Lazarian, & Vishniac 2002). 

Advances in computational power and numerical algorithms over the past few decades 
have made numerical simulations one of the important tools for quantitative investigations 
of magnetized turbulence (Shebalin et al. 1983; Cho et al. 2002; Cho & Lazarian 2003). 
However, turbulence usually covers a huge dynamical range from the macroscopic scales 
of turbulence generation to the microscopic dissipation scales. Current simulations have 
a dynamical range of a few hundreds to a few thousands and have not been able to give 
a complete description of energy flows in magnetized turbulence. These result in limited 
applications of them in astrophysics. Moreover, most of these studies arc limited to the 
MHD regime, where the background particles are strongly coupled with each other and 
can be treated together as a single fluid. The electron magnetohydro dynamics (EMHD) 
treats electrons and ions as two fluids. It is only applicable in a narrow frequency range 
between the electron and ion gyro-frequencies where the whistler dispersion relation is valid 
(Biskamp et al. 1999; Petrosian & Liu 2004). To address the heating of background particles 
by the turbulence, one has to assume that most of the background particles reach thermal 
distributions and arbitrarily extrapolate the turbulence spectrum into the dissipation range, 
where the MHD and EMHD formulism is usually invalid (e.g. Leamon et al. 1999; Tu et al. 
2002; Zhang & Li 2004; Wu & Yang 2006). Observations of solar flares, solar winds and 
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space plasmas, on the other hand, demand a detailed study of processes in the dissipation 
range, where background particles with different charge to mass ratios interact with the 
electro-magnetic fluctuations quite differently (Petrosian & Liu 2004; Liu et al. 2004, 2006; 
Lcamon et al. f998). This results in complicated wave dispersion relations, namely the 
dependence of the wave frequency on the wave vector (Andre f 985) and certainly affects the 
turbulence cascade. The couplings among waves and particles are even more complicated 
(Stix 1962; Xie 2004; Saito & Gary 2007). 

The diffusion approximation for the power spectrum in the wave number space has been 
a very powerful and efficient tool to study the turbulence cascade and dissipation over a large 
dynamical range. The 1-D models not only address the transition of turbulence from the 
large scale Kolmogorov phenomenology to the small scale IK cascade (Zhou & Matthaeus 
1990), but also are used to study the acceleration and heating of background particles by 
magnetized turbulence (Miller et al. 1995, 1996) and the damping of waves at small scales 
(Li et al. 2001). These studies have deepened our understanding of energy release processes 
during solar flares significantly. They, however, encounter difficulties in reproducing the 
broken-power spectrum of solar wind magnetic fluctuations (Li et al. 2001; Stawicki et al. 
2001). Due to the nonlinear nature of these 1-D diffusion models, the turbulence spectrum 
cuts off sharply at the wave number, where the cascade timescale becomes comparable to the 
wave damping time. The particle acceleration model also needs to be modified to reproduce 
the observed enhancement of high energy ^He ions during impulsive solar fiares since the 
1-D turbulence is damped by background ^He ions before it reaches the ^He gyro-frequency 
to accelerate the low energy ^He ions from a thermal background (Mason et al. 2002; Liu et 
al. 2004, 2006). 

Studies of weak Alfven turbulence with the quasi-linear treatment of electro-magnetic 
fiuctuations have shown that the cascade is anisotropic. It is also well-known that damping 
rates of different plasma wave modes by thermal background particles are very sensitive to 
the wave propagation directions. These anisotropics are critical to study energy dissipation 
through magnetized turbulence in a collisionless thermal plasma. A 2-D diffusion model is 
needed to advance of our understanding of magnetized turbulence and address the difficul- 
ties encountered with the 1-D models. Cranmer and Van Ballegooijen (2003) have showed 
recently that the heating of the background particles by Alfven turbulence are very sensi- 
tive to the 2-D and kinetic effects. To recover the critical balance and partially take into 
account the kinetic effects, they constructed a complicated diffusion-convection equation for 
the power spectrum with three dimensionless coefficients. A similar quasi-2-D model was 
proposed recently by Howes et al. (2007) to explain the broken power-law power spectrum 
of the solar wind turbulence. 
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In this paper, we study the general characteristics of the nonhnear 2-D diffusion model, 
the kinetic and damping effects of the Alfven-cyclotron fluctuations. In principle, all plasma 
wave modes, such as Alfven, fast and slow modes, should be included to have a complete 
description of the cascade and damping of magnetized turbulence. One then needs to solve 
a set of coupled diffusion equations for each wave branches (Andre 1985). This is a quite 
challenging task since the results will depend both on the interactions within each wave 
branch and couplings among different branches. Fortunately, Cho and Lazarian (2003) have 
shown that couplings among different wave branches are usually weak and these couplings 
decrease toward small spatial scales (Luo & Melrose 2006). Chandran (2005) showed that 
the couplings between Alfven and fast modes are significant only in the direction parallel 
to the large scale magnetic field, where the frequencies of the Alfven and fast modes are 
comparable. The couplings among different branches therefore can be separated from other 
processes near the dissipation range, where the kinetic and damping effects dominate. We use 
the exact dispersion relation for a cold plasma, which is a good approximation for the more 
general dispersion relation of a coUisionless thermal plasma. The diffusion coefficients are 
constructed as a function of the wave phase and/or group velocities so that the kinetic effects 
are treated self-consistently. The linear Vlasov equation is used to derive the thermal wave 
damping rates. Since the damping rate increases sharply with the increase of wave number. 
This quasi-linear treatment is expected to give a good approximation of the damping even 
for the nonlinear effect dominated quasi-two-dimentional component. The diffusion model 
actually does not distinguish the quasi-two-dimentional and wave-like fluctuations (Oughton 
et al. 2006). The balance between eddy turnover and wave propagation is revealed in the 
diffusion tensor. 

The Alfven-cyclotron branch has been studied extensively due to its simplicity and 
its prevalence in magnetized turbulence. It is chosen here to facilitate better comparisons 
of our model with previous results. In § [21 we discuss how the turbulence cascade can 
be studied using the diffusion approximation, which reduces the turbulence evolution to a 
nonlinear 2-D diffusion problem. The nonlinear diffusion equation can be solved numerically 
to obtain the power spectrum from the MHD region to the ^He gyro-frequency, where the 
Alfven dispersion surface cuts off. To compare with previous studies, we first presents the 
results for the cascade of Alfven turbulence § [31 where the dispersion relation for Alfven 
waves is adopted. The exact dispersion relation and damping rate are discussed in § [H The 
kinetic and damping effects are investigated in § [5], where we also consider its application to 
solar wind magnetic fluctuations. We discuss the implication of this model on the study of 
plasma heating and particle acceleration by magnetized turbulence, future work, and draw 
conclusion in § [Q 
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2. Kinetic Equation and Diffusion Tensor 

For weak or intermediate level turbulence, the most efficient way of energy cascading is 
through 3-wave interactions, with the resonance condition 

ki + k2 = k, u;(ki) +a;(k2) =u;(k), (1) 

where k and uo are the wave vector and the wave frequency, respectively, and c<j(k) gives the 
wave dispersion relation. For Alfven waves, a;(k) = v^kwi where is the Alfven velocity and 
k\\ is the parallel component of the wave vector. Since only oppositely directed wave packets 
can interact, this resonance condition further requires one of the interacting wave vectors, say 
k2, has to be perpendicular to the large scale magnetic, i.e., k2\\ =0 (Shebahn, Matthaeus & 
Montgomery 1983). As a result, Alfven wave turbulence cascades strongly in the direction 
perpendicular to the mean magnetic field. In the parallel direction, hydrodynamic eddy 
interactions prevent the wave cascading from vanishing but this can be a slow process. 

To study the wave cascade through the above 3-wave interactions, Goldreich & Shridhar 
(1995) first wrote down the integral form of the kinetic equation for the wave power spectrum 
>V(k,i) 



4 j d%d%(5(ki + k2-k)(ei-e2)''(k-e2)^>V(k2)[>V(ki)->V(k)]r(k2), (2) 



a>V(k, t) 
di 

where ej (j = 1,2) are the unit polarization vector of the specific Alfven waves. T(k) = 
r7(k)/[7y(k)^ + 4a;(k)^] is the timescale of wave-wave interaction determined by the "eddy 
damping rate" r] and the wave frequency. This equation strictly satisfies the 3-wave resonance 
condition and only allows energy cascade in the perpendicular direction. It does not take 
into account the nonlinear effects of turbulence and therefore cannot describe the general 
turbulence cascading in the 3-D wave vector space. By including Alfven-fast mode cross 
interaction, Chandran (2005) found that waves can cascade in all directions and solved the 
coupled kinetic equations for the Alfven and fast mode wave power spectra simultaneously. 
However, the Alfven wave spectrum does not appear to have any feature associated with 
"critical balance" first introduced by Goldreich & Shridhar (1994, 1995). 

In this paper, we focus on Alfven-Alfven wave interactions and extend the theory be- 
yond the MHD regime by using the general dispersion relation for the Alfven-cyclotron 
fluctuations. To incorporate the nonlinear effects caused by the eddy turnovers, we adopt 
the diffusion approximation, which leads to the general kinetic equation: 



A.|-W(k,t) 



r(k)>v(k,t)-^i|, (3) 



8 



where the terms on the right-hand-side represent the wave generation, cascade, damping and 
leakage processes, respectively. In the following, we will ignore the leakage term and treat the 
source term as a 5-function at certain large scale and study how the wave spectrum depends 
on the cascading and damping processes. By proper choices of the diffusion tensor, the 
cascade caused by the wave-wave resonances and eddy turnovers may be modeled (Zhou & 
Matthaeus 1990). We are interested in the solution that decreases rapidly with the increase 
of the wave number. Equation ([2j) then shows that contributions of the 3-wave interactions to 
the evolution of the power spectrum is dominated by processes with k2 1. This justifies the 
diffusion approximation for the 3-wave interactions. For k2 ^ 1, (Wi — Wk) — k2 • VVV(k), 
which corresponds to a flux in the wave vector space. Comparing equations ([2]) and ([3]), it 
is obvious that the amplitude of the diffusion tensor is proportional to the total turbulence 
energy. 

Following Zhou & Matthaeus's (1990) theory, we construct the diffusion tensor with the 
nonlinear timescale tnl associated with eddy turnovers and wave crossing time ta. The eddy 
turnover time can be estimated with tnl = ivk)~^, where the eddy velocity v ~ (Wk^Y^'^. 
For the wave crossing time, one must consider the anisotropy of Alfven-cyclotron dispersion. 
Near the perpendicular direction, since ^2,11 — 0, any wave can always interact with this 
large scale wave, i.e. ta ^ 00. The wave crossing effect can be ignored, and the cascade is 
dominated by the eddy turnover; In the other directions, wave packets with a size of ~ l/fc|| 
cross each other at the Alfven speed. In general, we have ta — (fAfc||)^^Ill Then the 3-wave 
coupling time = (t^l +'%^)~\ and the wave cascading rate r^J. ^ 't^I'Ts- Then the locally 
isotropic diffusion tensor 



where 6ij is Kronecker's delta with i,j indicating the three bases of the 3-D wave vector 
space, C is a dimentionless cascading constant corresponding to the Kolmogoroff constant 
for hydrodynamic turbulence, and 6 is the angle between the wave vector k and the mean 
magnetic field. 



In the strong turbulence limit, v ^ va, this diffusion tensor recovers the isotropic 
Kolmogoroff cascade. However, in the weak turbulence limit, i.e. v <^ va, it does not 



^Note that, following Kraichnan's (1965) argument for an isotropic Alfven wave turbulence, Zhou & 
Matthaeus (1990) obtained r4 = (vAk)^^ for their 1-D diffusion model. 




(4) 



3. 



Cascade of Alfven Turbulence 
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reproduce the isotropic Kraichnan spectrum due to the global anisotropy of the diffusion 
coefficient caused by the anisotropic Alfven wave crossing time. Since the Alfven wave 
crossing time ta is much longer in the perpendicular direction than in the parallel direction, 
the wave cascading rate is large in perpendicular direction and decreases dramatically with 
the increase of k\\. The wave energy contours in Figure [1] show the numerical result of such 
anisotropic cascading. The turbulence energy cuts off at the cascading front. In fact, if we 
assume a Kolmogoroff power spectrum, the Alfven wave crossing process starts to dominate 
and suppress the wave cascading aXk\\ > k^l^ , which is similar to the critical balance relation 
introduced by Goldreich & Shridhar and explains the k\\ ~ scaling of the cascading front. 
This cascading front in wave vector space extends until it reaches the numerical or physical 
boundary (damping or non-MHD effect at large k) . However, although the cascading energy 
flux is highly concentrated on the perpendicular direction at small k\\, the diffusive process 
still carries wave energy to large k\y By assuming an infinite damping at certain large k 
and let turbulence evolve freely at any k below this infinite damping boundary, we can see 
(Figured]) the turbulence fills all possible wave vector space in our numerical simulation. In 
our model, the energy transfer to large through two mechanisms: the reverse cascading 
on perpendicular direction and the possible slow cascading in parallel direction. For MHD 
turbulence, many publications () consider the cascading in parallel direction void beyond 
certain critical value, i.e. k\\ ~ k'^/_'^ . However, there is no theory against inverse cascading 
on perpendicular direction. As a result, we argue that our diffusion theory is valid and the 
wave energy fills all possible wave vectors before it gets damped. 

The "steady state" spectrum we obtained based on this simple diffusion model is a 
quasi-isotropic Kolmogoroff- like spectrum, with highly anisotropic energy flux concentrated 
on perpendicular direction. This is due to an "isotropic" diffusion tensor (see Figure [T]). 
Theoretically, by setting the wave cascading rate to be large at small k\\ and vice versa, 
we are able to simulate anisotropic turbulence energy flux concentrated on perpendicular 
direction. However, since we set the diffusion coefficient to be the same on parallel and 
perpendicular direction (-D||,|| = D^ i_) for at any k, we can only reach an isotropic spec- 
trum given enough turbulence evolving time. An apparent improvement of this model is to 
construct the parallel and perpendicular terms of diffusion tensor differently. We can take 
parallel and perpendicular direction as the diffusion tensor's principle axes (eigenvectors) 
and the corresponding diffusion coefficients as its eigenvalues. For cascading on perpendic- 
ular direction, since A;2,|| = 0, the wave cascading is limited by wave eddy turn over time 
only, then we can write the -D±,± = k'^T'^^g ^ = k^r^]^. For cascading on parallel direction, 
small scale wave ci;(ki) collides with large scale wave uj{]<.2) and gradually obtain a drift in 
wavelength (ki —>■ k). Although this drift can be both way (to smaller or larger scales), 
the energy gradient in wave vector space keeps the energy flow from large to small scales 
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stronger than the inverse. Since we are interested in the change of wavelength during a sin- 
gle wave period of ki, its own frequency a;(ki) limits the cascading rate. Thus we combine 
this wave frequency with eddy turn over rate to obtain cascading rate on parallel direc- 
tion T^~^|| = t^\/{t^\ + Ta^)- And the corresponding diffusion tensor can be written as 
Z^ll^ll = k'^T^Jg II • Thus, we obtain the form of the anisotropic diffusion tensor as 

where t^]^ and r^^ is defined the same way as in isotropic diffusion model. The resulting 
spectra in different directions (see Figure [2]) show such an anisotropy. Note that for both 
isotropic and non-isotropic diffusion tensor model, we have larger cascading rate at small 
and the energy flux is highly concentrated on perpendicular direction. 

Instead of exploring more possibilities on the forms of diffusion tensor, we refer to 
observation of turbulence spectrum in solar wind to point us the direction. The discussion 
can be found in § [51 



4. Dispersion Surface beyond MHD Region and Thermal Damping 

The Solar Wind observation suggests that the Alfven wave turbulence cascade to smaller 
scales than the MHD regime. Along parallel direction (for all kj_ < 1/pp, where Pp is proton 
Lamour radius), the Alfven wave dispersion surface diverge from simple u ~ kvA cos 6 around 
kvACOsO < 0.5Qp and flatten into He cyclotron oscillation {u = 0.5Qp) for larger k\\s. We 
discuss how this affects the wave energy cascading below. On the perpendicular direction, 
the dispersion surface bends into kinetic Alfven wave (KAW) at kpp ~ 1. The effects of 
KAW dispersion relation and its high damping rate on turbulence cascading is even more 
complex. However, as we will discuss in § [5], in typical solar wind conditions (where plasma 
beta /3p < 1 or ~ 1), KAW is Doppler-shifted to a high frequency around or beyond the 
satellite observation limit. Thereafter, we limit our study to a wave vector space below 
KAW range in this paper. 

When the damping rate is small, the hot plasma dispersion surface can be calculated 
through numerically solving Vlasov equation, where Lauren series is used to approximate 
the Z-function in the solution. Because of this approximation, when the damping rate is 
getting close to the wave frequency from below, this numerical method fails and there is no 
simple scheme that can guarantee a solution. Figure H] shows the numerical result of Alfven 
wave dispersion surface using Waves in Homogeneous, Anisotropic, Multicomponent Plasmas 
(WHAMP) code (Ronnmark 1982). The missing segments in the dispersion relation are 
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the wave vectors where the numerical solution fails due to the relative large damping rate. 
Fortunately, as we can also see from Figure 5, the cold plasma approximation provides a 
close fit to the "real" dispersion relation in the wave vector range of our interest (i.e. beyond 
MHD but below KAW). Since we have the analytical form of cold plasma dispersion relation 
(see Appendix |A]) , we can solve the smooth function wave frequency co'(k) and its gradient 
(i.e. wave group velocity) Vg = Vki^(k) at all the wave vector grid points for the use of 
our numerical simulation on turbulence energy cascading. Thereafter, we use cold plasma 
dispersion relation to study this diffusion process in the following. 

It requires a detailed study to understand how the dispersion surface beyond the MHD 
regime affects wave cascading. In the region where the dispersion surface starts to diverge 
from MHD approximation, i.e. k\\VA ~ 0.5Qp, this surface, like in MHD regime, is flat 
on perpendicular direction (i.e. duj/dk± <^ duj/dk\\). It's easy to demonstrate that the 3- 
wave resonance condition (Equation ([1])) still requires the smaller one of the interacting wave 
vectors ^^2 to be in perpendicular or near perpendicular direction. So, at region k\\VA ~ 0.5f2p, 
the Alfven wave cascading is still highly preferred in perpendicular direction and we can 
simply adapt the wave crossing time as = Vg ■ k ~ duj/dk\\k in our diffusion tensor. In 
the region far beyond MHD, i.e. k\\VA ^ O.Sfip, the cold plasma dispersion surface flattens 
out in all directions and the Alfven wave packet becomes stationary. It is not obvious how 
the turbulence will cascade in this region. In one aspect, three wave resonance condition still 
requires one wave vector to be perpendicular, i.e. we need the cascading rate (or diffusion 
tensor) to be larger in perpendicular direction. However, Vg ■ k = can neither give us such 
a parallel cascading rate nor the wave crossing time, instead = l/u;(k) provides a better 
approximation. In the other aspect, turbulence cascading can be as simple hydrodynamic 
when the oscillation is stationary, with eddy turn over rate serve as the only limit on wave 
cascading. Then = Vg ■ k = is just what we need for the diffusion tensor. Thereafter 
we try two different approaches to approximate the wave crossing time: ta = (vg ■ k)"-*- or 
ta = uj~^. The final spectra based these two different wave crossing time are shown in [3J 
They both contain a spectral break at k\\VA ~ ^a-, where the wave dispersion surface deviats 
from the MHD approximation. These breaks, though appealing, only have an index change 
less than one. Thereafter neither of these breaks can explain the observed broken power-law 
spectrum in solar wind for Alfven turbulence (Leamon et al. 1989). Furthermore, as we 
show in § [5l these two wave crossing rate models can can hardly produce any observational 
difference in in solar wind when thermal damping is taken into consideration. Thereafter, 
we leave this theoretical work to further study and use ta = in our simulation. 

To summarize, we get the simple isotropic form of diffusion tensor Dij for MHD and 
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beyond: 



- X ^ p ^NL - ^ ^ (a) 



The form of anisotropic diffusion tensor remains the same as Equation ([5]), only that 



The wave damping rate, can only be obtained through solving Vlasov equation. For 
parallel and quasi-parallel propagating waves, Swanson (1982) simplified the Dielectric Ten- 
sor and obtained the damping rate for weak damping approximation by electron cyclotron 
oscillations. To study all the damping effects we generalize the formular to calculate the 
cyclotron damping from all species of particles: 



UJi 



^ \^^p,s r f 

y -^exp - 



kv. 



2 

+ > ^-T^^v 



LUr 2k sr^ 2uj — t^Vtg 



where different s stands for different particles. In our study, s ranges from 1 to 3 standing 
for electron, proton and ^He respectively. Note that although this formula gives the value 
of UJi all or k, the approximation fails at cjj ~ ^c,s- At these region when approximation 
fails, Swanson (1982) shows that cjj ^ Ur- This result is good enough for our study of wave 
cascading and damping, because the exact value of damping rate is not needed for strong 
damping. 

For wave propagating in other angles than parallel, we use WHAMP code to solve the 
damping rate numerically. Similarly, when tu, ~ tu^, the weak turbulence approximation is 
not valid and the code fails to converge. At this region, we roughly estimate the damping rate 
with power-law extrapolation (linear in our log-log graph) since the exact value of damping 
is not necessary when it is strong enough to cut off the wave. 



5. Numerical Results and Observation 

We solve the 2-D time dependent wave kinetic equation with Alternative Direction 
Implicit (ADI) scheme on a log-log uniform grid to study this wave cascading process. A 
reflective boundary condition is used at large scale boundary fcmm- For small scale bound- 
ary, we set kmax large enough that all wave energy is damped to zero. We put a small 
constant injection at large scale starting from t = to simulate the turbulence excitation, 
i.e. Qy\)(k,t) = QoS{k — ko)6'(t) in wave kinetic equation (Equation ([3])). 



5.1. Anisotropy by Isotropic Diffusion tensor 



The numerical result of steady state turbulence spectrum for isotropic diffusion tensor 
is highly anisotropic (see Figure [5]). To understand the result from cascading and damping, 
we perform simulation without thermal damping (instead of including thermal damping at 
proper k, we set an artificial isotropic cut-off at very large k). In such a case, the injected 
turbulence evolves into a quasi-isotropic Kolmogoroff spectrum in MHD region and breaks 
slightly beyond that. The wave cascading rate, as we discuss on MHD turbulence in § [21 is 
large near perpendicular direction and small when otherwise (see Figure [5] and Figure [6]). 
Meanwhile, the wave damping is weaker near perpendicular direction than in the parallel 
direction. This anisotropic damping, when combined with anisotropic cascading, makes the 
Alfven wave spectra cut off at different magnitude of k on different angles, and the cut-off 
region extends over one order of magnitude in wave vector space when integrated over over 
the fc-sphere. As shown in Figure El the extended cut-off region has a spectrum close to 
a power law with index 3 ~ 4. This spectrum suggests us to use it in explanation of the 
observed broken power law spectrum in solar wind Alfvenic turbulence. To compare with 
the observation we integrate the steady state turbulence over wave vector space to obtain 
the spectrum in Doppler-shifted frequency space of the spacecraft (Leamon et al. 1999): 



where / is the spacecraft-frame frequency, F{f) is observed spectrum and is the Dirac 
delta function. As shown in Figure [HI after this Doppler-shifting integral, the extended 
cut-off region becomes a power-law-alike spectrum. And this spectrum provides a good fit 
to the spectrum of the Alfvenic oscillation in solar wind observed by Leamon et al. 1998. 
Furthermore, since this broken power law is a direct result of the extended cut-off region in 
wave vector space, we can approximate the break frequency with the starting point of 
this cut-off region. If we define the anisotropic cut-off points in wave vector space as kc(0), 
then the observed break frequency can be roughly estimated by 



We can see from Figure [6] that kc is much large at near perpendicular direction, thus we 
can expect that kc(6') ■ Vgw increase dramatically with the increase of Qbv- However, the 
ming{...) in the formula mixed up the simple relationship and we obtain a weakly dependence 
of the Ubf on Qbv (see Figure [TUI) . This result qualitatively agrees with two similar cases 
observed by Leamon et al. (1998), where all parameters are close except for the angle of 
solar wind. It may come to attention that in the Qbv = 87° case, the observation, the 
second index of the broken power law appears much harder than the Qbv = 23° case, which 





Ubf min(kc(6') ■ Vgw) 



min(/cc,||V5H/ cos{Qbv) + kc,±Vsw sin(eBy)) 
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is on the contrary to our numerical result. We argue that in the Qbv = 87° case, the second 
power law segment is too close to the background noise to fit an accurate index number. 
Thereafter, we only need to take the break frequency into serious consideration. We further 
discuss this effects of noise saturation on fitted index number in § 15.31 

We emphasize that this broken power law spectrum requires anisotropy of cascading 
as well as damping. Figure [TO] gives a direct comparison between the kc(6') and equal Ui 
contour in wave vector space. Note that the Ui contour would mark the cut-off in different 
direction if the cascading was isotropic. Since Ui is more isotropic, the cut-off region would 
not be extended enough to make a broken power law spectrum without anisotropic cascading. 
The numerical test with an artificial isotropic cascading (we use Kolmogoroff cascading with 
lower injection rate) confirms our conclusion (see Figure [7]). B 



5.2. More Anisotropy by Anisotropic Diffusion tensor 

The non-isotropic diffusion tensor provides a Kolmogoroff index spectrum in perpendic- 
ular direction only. In other directions, the energy spectrum is a much steeper. Although 
the spectral cut-off spreads to a even larger range of k in such a case, the energy in parallel 
direction contributes very little to the /c-sphere integrated spectrum. As shown in Figure [HI 
there is no extended cut-off region for this model. As a result, when integrated into Doppler- 
shifted frequency space, the spectrum after break is much more steeper than observation 
range (with index 7) even we fit it as a power law. This disagreement with observation is 
also found in many numerical or analytical studies (Li 19XX, Howes et al. 2007). § In a 
word, a simple thermal damping can not produce the observed broken power law in a 1-D 
cascading model (or 2-D if we take the two orthogonal direction in perpendicular direction 
separately). When most of the turbulence energy is on perpendicular or close to perpen- 
dicular direction, we can only refer to change of cascading rate itself to produce a broken 
power law. The KAW, with a dispersion relation diverges from the cold plasma or MHD 
approximation, provides a different cascading rate. However, Howes et al. (2007) show that 



^By assuming a hard spectral index (72 = 5 in the paper) beyond certain tOi contour in wave vector space, 
Leamon et al. (1999) obtained a broken power law spectra in frequency space. They argue that thermal 
damping alone can explain the observation. However, our simulation shows that the turbulence spectra cut 
off very steeply when damping rate is large enough. Thereafter their assumption is invalid, so does the 
conclusion. 

■^Howes et al (2007) suggest that the spectral cut-off, when combined with a observational noise saturation, 
may appear to be a power law spectrum. In this paper we stick on the broken power law conclusion observed 
by Leamon et al. (1998) and make a direct comparison between observation and our theory. 
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the spectral index of KAW is only -7/3, much less than second index of the observed broken 
power law. Furthermore, KAW only starts to appear at /c_LPp ~ 1 {pp is proton Larmor 
radius), which is Doppler-shifted to a spacecraft frequency 

~ kperp ■ Vsw ~ ik_i_pp ^^^'^ p-^^^)up (9) 

where Up is proton gyrofrequency. In typical Solar Wind conditions, where /3p < 1 or ~ 1 and 
Vsw ~ lOf A, this frequency is 10 times the proton gyrofrequency or over. Whereas most of 
the observations (Leamon et al. 1998 & 1999, Bale et al. 2005) contains a break frequency 
smaller than lOup. Although the available observation results are not complete enough to 
rule out KAW as the reason for the broken power law, a different mechanism other than 
KAW is suggested. 

To generalize, we are unable to reproduce the observed broken power law spectrum with 
our anisotropic diffusion tensor model. We also show that this difficulty persists on any model 
that producing a highly anisotropic turbulence spectrum. The isotropic diffusion tensor 
(with non-isotropic cascading rate), on the contrary, provides a good fit to the observation. A 
further observational evidence is the weak angle dependence of the break frequency (Leamon 
et al. 1998), which can be very difficult to explain for a highly anisotropic spectrum but can 
easily fit into the isotropic diffusion tensor model. Thereafter, in the following discussion 
and observational prediction, we use the isotropic diffusion tensor model only. 



5.3. Break frequency and observables 

In our model, the broken power law is produced by the wave spectral cut-off at different 
magnitude of k in different direction, hence the breakpoint in spacecraft-frame frequency 
space corresponds to the cut-off in wave vector space indicated by Formula ([8]). This model 
allows us to draw some conclusions on the relation between spectral breakpoint and the 
observable parameters. 

We categorize the observables into four independent parameters, Alfven Mach number 
Ma, plasma beta Pp, Vsw, and Qbv- Obviously, a faster solar wind will give a higher 
break frequency. This relation is defined by the Doppler-shift term k • Vsw in Equation 
([7]) rather than depends on turbulence model. Instead, we are interested in relations that 
are turbulence model dependent as below. First, since the cut-off (or break in observation) 
is produced by turbulence damping overwhelming its cascading, the relative value between 
turbulence intensity and thermal damping rate will affect the breakpoint the most. We can 
expect that a high Alfven Mach number will give a high break frequency and high plasma 
beta will reduce break frequency. Our numerical simulation in Figure [11] confirms these 
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relations. Note that since the cascading rate still contains one undetermined free parameter, 
the cascading constant C in Equation (jl]), we can only predict the z/b/-M^ or Uhf-jSp relation 
rather than provide the exact value of break frequency from observed turbulence parameters. 
Although it is interesting to notice that C = 1 gives a very close match to the observation 
made by Leamon et. al. (1999) in both break frequency and turbulence intensity, we @ 
Secondly, the projection effect from Qbv also affects the break frequency. It is obvious that 
the break frequency increase at larger projection angles. However, the angle dependence is 
not as critical as predicted by any 2-D (wave energy cascades only to perpendicular direction) 
turbulence model. With 3-D wave spectrum and cut-off wave vector kc in our model with 
isotropic diffusion tensor, the minimum value of kc • Vsw depends on both the shape of 
kc contour and the angle of solar wind. A more isotropic cut-off wave vector will give less 
dependency on the angle of solar wind. On the other hand, if most of the wave energy resides 
in perpendicular direction, as suggested by critical balanced spectra, the kc contour will be 
highly prolonged in perpendicular direction and the angle of solar wind will strongly affect 
the break frequency. The observation by Leamon et al (1998, Figure 6) fits better with a 
more isotropic model as ours but further observation is required to rule out any one. 

The spectral index after the break, on the other hand, is also related to these observable 
parameters as we show in Figure [10] and [TTl However, since the broken power law is only 
an observational effects of the extended cut-off region rather than a mathematical secondary 
power law spectrum, the fitted index number highly depends on the low energy cut-off of the 
observation, and the noise saturation also effects the result (Howes, et. al. 2007). Thereafter, 
we do not attempt to compare theory with observations unless further detector can measure 
a definite result on the index number. 



6. Conclusion 

With diffusion approximation, we are able to construct a diffusion tensor discribing 
Alfven-cycloton turbulence cascade. This cascade model, when coupled with thermal damp- 
ing, provides a broken power law spectrum after integration, which matches the observed 
solar wind plasma turbulence spectrum by Leamon et. al. (1999). By comparing the 
model with observation, we conclude that the observed broken power law spectrum comes 
from a quasi-isotropic turbulence spectrum that cuts off at different k in different direc- 



^In principle, we are able to determine the cascading constant by fitting our theoretical breakpoint with 
this observation. However, we are not able to obtain all the observational parameters for the individual case 
from Leamon et. al. (1999). We pick a typical Alfven wave velocity in solar wind as va = 40km/s, and use 
this value for all our simulations. 
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tion. Based on this model, we predict the observed break frequency is propotional to Alfven 
Mach number, and antipropotional to plasma beta. This result can be subjected to further 
observational tests. 

This work was carried out under the auspices of the National Nuclear Security Ad- 
ministration of the U.S. Department of Energy at Los Alamos National Laboratory under 
Contract No. DE-AC52-06NA25396. This research was partially supported by NSF grant 
ATM-0312344, NASA grants NAG5-12111, NAGS 11918-1 (at Stanford). 

APPENDICES 



A. The Exact Solution of Cold Plasma Dispersion Surface 

The equation for cold Plasma dispersion surface can be constructed through solving 
Maxwell's equation for plane waves (Stix 1962). A nontrivial wave solution to the equation 
requires the refractive index n = ^ to satisfy equation, one obtain the wave equation, 

nx(nxE) + K- E = (Al) 

where K is the dielectric tensor defined by 

J - iojeoE = -iLoeK ■ E. (A2) 

By defining P, R,L,S ^ ^f^iR + L),D^ ^f^iR - L) as & function of lo, and A, B, C terms 

as a function of P, R, L and 9, Stix (1962) was able to write all the wave solution in the 

B ± F 

form ri^ — — ^-^ — . Although this solution is segmented and has may poles at particles' 
cyclotron frequency, each physical mode of dispersion surface is both continuous and smooth 
almost everywhere. In the following, we express these continuous modes by combining with 
segmented solutions of the equation (??), i.e. 

Alfven branch. 



k = ±uj^^^ u;e[0,Qa) (A3) 

R and L term reach their first pole {fla), and we have lim^_^Q- k — oo. This point is the 
end of Alfven branch. 

Fast branch, 

k = 



W^^2J^ ^e[o,Q-] 
1 jp 
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Similarly, Fast branch ends at Qp with /c — > oo. Although is a pole for R and L, 
lim^^j^- = lim^^j-2+ ^ oo. To understand the switching sign at Qa, I study the 
simplified problem at ^ = 0. In such a case F — 2\J P'^D'^ — 2 | PD \ and I get 

/ I PD I \ 

k ^ ±uj \S ±^-^^\ ^±u;{S±sign{PD)Df^^ (A4) 



Since P < for all uj G [0, 1], and D switch sigh at 0„, the discontinuity is only introduced 

B±F 
2A 



by our attempt to write k{uj) in an explicit form. By switching sign in at D = 0, I can 



follow the continuous dispersion surface. 
Whistler branch 



k 



I B + F 



' B — F 



2A 

Whistler branck starts at uji, which solves the equation R{uj) =0. ci;i ~ 0.585fip (very 
unsensitive to density and field strength) by numerically solving the equation. P = at Ue, 
the fomular reach the pole in parallel direction = 0), and the Whistler branch ends at the 
electron Langmuir oscillation. 
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Fig. 1. — Wave energy W{k\\,k±) contours at intermediate (left panels) and steady state 
(right panels) for weak turbulence without damping. The upper panels show the result 
of isotropic model (-D||,|| = D±^±) and lower panels show the result of anisotropic model 
7^ D±^^. The red line labels the critical balance relation introduced by Sridhar & 
Goldreich (1994). 
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Fig. 2. — Power spectra of Alfven wave turbulence for two diffusion tensor models. The 
result from isotropic model (-D|| || = -D±,_l) is shifted by 10 times. The isotropic model 
provides a close to isotropic Kolmogoroff spectrum and the anisotropic model D\\^\\ ^ D±^± 
provides a steep spectra in directions other than perpendicular direction. 
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Fig. 3. — Wave energy spectra with two different wave crossing rates beyond MHD regime. 
The breaks at k\\VA ~ 0.5Qp in both spectra are due to the dispersion relation co'(k) devi- 
ating from the simple MHD relation uj = k\\VA- The two spectra, as discussed later, can 
not be differentiated by observation because of the strong thermal damping for solar wind 
conditions. 
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Fig. 4. — Plasma dispersion relations under typical solar wind conditions comparing the 
one under cold plasma approximation. The former one is solved numerically with WHAMP 
code. 
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Fig. 5. — Left: Wave cascading rate (solid) vs. damping rate (dotted) in different direction 
(with different colors). In any direction, the crossing point of two rates gives a rough es- 
timation of where the energy spectrum cuts off. Right: The resulting spectral cut-off that 
spreads into 1.5 order of magnitude in different direction, which can produce a broken power 
law spectrum when integrated over 6. 
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Fig. 6. — Left: Green and red lines show the sample contour of cascading rate and damping 
rate. When the two contours with equal value cross (the two contours cross at 45° in 
this example), we obtain the critical kc, which roughly labels turbulence spectral cut-off in 
different directions. The cascading rate increase at large angles and vice versa, so the kc 
goes above the damping rate contour at large angles and below the contour at small angles. 
This variation of k^. over angle produces an extended cut-off region in k space when the 
turbulence spectrum is integrated over angle (see Figure [8]). Left: The resulting k^ curve for 
different conditions. Solid, dash and dotted lines show the case where (M4 = 0.6, (3p = 0.41), 
(M4 = 1.1, Pp = 0.41) and (M4 = 0.6, /3p = 1.0) respectively. Apparently, the kc contour is 
shifted to high ks at high turbulence intensity and lower ks at high damping rate. 




Fig. 7. — Spectral cut-off in different directions. A Kolmogoroff (Homogeneous) cascading 
with anisotropic damping is simulated (the middle bunch of lines) to compare with result 
our diffusion tensor model. 
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Fig. 8. — Panel (a): Wave turbulence spectrum as function of k after integration over angle. 
The result from isotropic diffusion tensor model (red) steepens slowly and the cut-off region 
extends about one order of magnitude in wave vector space; the result from anisotropic 
tensor model (blue) cuts off rapidly around 1/pp {vA f=i Vth,p in this simulation). The 
latter spectrum agrees with Howes et al. 2007. Panel (h): Wave turbulence spectrum as 
function of wave frequency when solar wind Doppler-shift is absent. Both spectra show a 
Kolomogoroff index and cut off at Hellium gyrofrequency (i.e. the upper limit of cold plasma 
Alfven dispersion surface). Panel (c): The spectra in Doppler shifted space-craft frequency 
for two Dij models. The result of isotropic tensor model can be fit into a broken power law 
spectrum, with 71 = —1.67, 72 = —3.1 and break frequency Ubf = 0.2Hz. 
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Fig. 9. — Observed interplanetary turbulence power spectrum (solid) vs. simulation results 
(dash). The data are extracted from Figure 1. by Leamon et. al. (1999). For this observation 
period, B — 6.3nT, Pp — 0.48, Qbv — 38°, Vsw — 517km/s, va is unknown (we use 40km/s 
for our simulations). With these parameters as input, our isotropic model (-D||,|| = -D±,±) 
provides an excellent fit to the observation. Leamon et. al. (1999) fit the observation 
with a broken power law with indices 71 = —1.67, 72 = —2.91 and break frequency at 
Ubf = 0.235Hz. Our simulated spectrum, when fit to a broken power law down to 10~^PSD, 
provides 71 = —1.67, 72 = —2.97 and break frequency u^f = 0.200Hz.) 
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Fig. 10. — Numerical simulation result for angle dependence of break frequency. Left: The 
numerical simulation of spacecraft frequency spectra at different solar wind to mean magnetic 
field angle, Qbv, with all other parameters equal (Ma = 0.6 and Pp = 0.41, see Figure [9] for 
values of specific observables) . Broken power law fits provides break frequency, Uhj = 0.15 
for Qbv = 15° and u^j = 0.24 for Qbv = 60°. Right Setting observational turbulence 
power cut-off as PSDc = 5 x 10-^(nT2Hz-^) (blue) or PSDc = 1 x 10-^(nT2Hz-^) (red), 
we simulate the relation between break frequency and solar wind angle in the upper panel 
and relation between after-break-index and solar wind angle in the lower panel. (Note that 
when observational noise and saturation is taken into account, the index after break may 
not be accurately measured.) Numerical simulation shows that the break frequency weakly 
depends on the solar wind angle. The other parameters are fixed at the same value as in 
Left panel. 
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Fig. 11. — Left Same as Figure fTOl but dependence on the Alfven Mach number. Other 
parameters are fixed at the same value as in Figure [H i.e. Qbv = 38°, Pp = 0.41. Right: 
Dependence on plasma beta. Other parameters are fixed at the same value as in Figure [H 
i.e. Bbv = 38°, Ma = 0.6. 



